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Abstract. We present a recently developed code for detailed NLTE calculations of accretion disc spectra of cataclysmic 
variables and compact X-ray binaries. Assuming a radial structure of a standard a-disc, the disc is divided into concentric rings. 
For each disc ring the solution of the radiation transfer equation and the structure equations, comprising the hydrostatic and 
radiative equilibrium, the population of the atomic levels as well as charge and particle conservation, is done self-consistently. 
Metal-line blanketing and irradiation by the central object are taken into account. As a first application, we show the influence 
of dilferent disc parameters on the disc spectrum for the helium cataclysmic variable AM CVn. 
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1. Introduction 

Accretion discs are components of objects as diverse as proto- 
planetary systems, active galactic nuclei, cataclysmic variables 
or X-ray binaries. A high fraction of the luminosity of these 
systems may be generated by the accretion disc itself. To under- 
stand these objects and interpret the observational data a model 
' of the accretion disc as reliable as possible is necessary. The 
. aim of our work was the development of a program package 
for the calculation of synthetic spectra and vertical structures 
of accretion discs considering the physical processes in the disc 
as accurately as possible. 

A fully three-dimensional radiation-hydrodynamic treat- 
ment is presently still impossible because of the enormous 
numerical costs. In the case of a geometrically thin c-disc 
(Shakura & Sunyaev 1973), where the disc thickness is sig- 
nificantly smaller than the disc diameter, the radial and verti- 
cal structures can be decoupled. Under the assumption of ax- 
ial symmetry and by dividing the disc into concentric rings 
the determination of the vertical structure becomes a one- 
dimensional problem. The dissipated energy in each disc ring 
is radiated away at the disc surface, the energy flux can be ex- 
pressed as effective temperature. As an approximation, it can 
be assumed that the disc rings are radiating like black bodies. 
An improvement of the models can be obtained by describ- 
ing the disc rings by stellar atmosphere models of the same 
effective temperature, see e.g. Kiplinger (1979), Mayo et al. 
(1980) or La Dous (1989) in the case of cataclysmic variables 
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and Kolykhalov & Sunyaev (1984) or Sun & Malkan (1989) in 
the case of AGN. Unfortunately, neither black bodies nor stel- 
lar atmosphere models reproduce spectra of accretion discs in 
an adequate manner (Wade 1988). Meyer & Meyer-Hofmeister 
(1982), Cannizzo & Wheeler (1984) and Cannizzo & Cameron 
(1988) calculated the radiative transfer using the diffusion as- 
sumption. This assumption is only valid at large optical depths, 
but the spectrum is generated at optical depths around r ~ 1, 
where neither the diffusion assumption nor the assumption of 
local thermodynamic equilibrium (LTE) is fulfilled. Only solv- 
ing the radiative transfer equation self-consistently with the 
structure equations allows the calculation of realistic accretion 
disc spectra (Kriz & Hubeny 1986, Shaviv & Wehrse 1986, 
Stoerzer et al. 1994). In the last decade much work in this field 
was done e.g. by Hubeny & Hubeny (1997, 1998) and Hubeny 
et al. (2000, 2001). 

Following the path mentioned above, we have developed 
our program package AcDc (Accretion Disc code) for the de- 
tailed calculation of vertical structures and NLTE spectra of 
accretion discs. For each disc ring the equations of radiative 
and hydrostatic equilibrium as well as the NLTE rate equa- 
tions for the population numbers of the atomic levels are solved 
consistently with the radiation transfer equation under the con- 
straint of particle number and charge conservation. Full metal- 
line blanketing as well as irradiation of the accretion disc by 
the central object can be considered. By integrating the spectra 
of the individual disc rings, one obtains a complete disc spec- 
trum for different inclination angles, where the spectral lines 
are Doppler shifted according to the radial component of the 
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Kepler rotation. This is shown in Sect. 2, whereas in Sect. 3 
we show first applications of the developed program package 
to examine the influence of different parameters on the vertical 
structure and the spectrum of an accretion disc model for the 
helium cataclysmic variable AM CVn. 

2. Vertical Structure of Accretion Discs 

We assume a geometrically thin (disc thickness is significantly 
smaller than the disc diameter) stationary accretion disc (a- 
disc, Shakura & Sunyaev 1973). We also assume that the mass 
of the disc is much smaller than the mass of the central object, 
so we can neglect self-gravitation. Introducing the surface mass 
density E as 



side view 



H/2 

- 2 J pdz. 



(1) 



with mass density p, geometrical height z above the midplane 
and total disc height H, the radial dependence of E reads fol- 
lowing Shakura & Sunyaev (1973) 



(2) 



Here, R denotes the distance from the central object, ^* the 
radius of the central object, M the mass accretion rate and v the 
kinematic viscosity, defined as 



V - aCsH , 



(3) 



with sound speed Cs and the parameter a being a measure of 
the efficiency of angular momentum transport through the disc. 
The radial distribution of the effective temperature T^ff then, 
following Shakura & Sunyaev (1973), can be described by 



3GM^M 



-(r 



1/4 



(4) 



with denoting the mass of the central object, G the gravita- 
tional constant and cr the Stefan-Boltzmann constant. The ac- 
cretion disc is divided into a set of concentric disc rings (cf. Fig. 
[TJ. For each ring the vertical structure is calculated by solving 
the set of equations described in the following two subsections, 
assuming a plane-parallel geometry. 

2.1. Radiation Transfer Equation 

In order to compute the radiation field, which determines the 
atomic population numbers, the radiation transfer equation has 
to be solved. This equation describes the modification of the 
specific intensity ly of a ray due to absorption or emission along 
its path ds through the accretion disc (cf. Fig.|5J. The radiation 
transfer equation for the shown geometry then reads 

dly(v,fi,z) 



dz 



-Xviv, z) Iy(v, fi, z) + r]y{v,z). 



(5) 



with the absorption coefficient Xv and the emission coefficient 
rjy. Introducing the source function Sy 



Sy^'-^ 

Xv 



(6) 
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Fig. 1. Geometry of the accretion disc, divided into concentric 
rings. 
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Fig. 2. Radiation transfer through the accretion disc layers. 

the radiation transfer equation is 
dly{v,n,z) 



dz 



■Xy{v,z){Iy{v,n,z) - S y{v, fi, z)) ■ 



(7) 



The solution of this equation is a complicated problem, because 
the source function depends on the radiation field itself. Within 
an iteration scheme (see below), one solves the radiation trans- 
fer equation formally assuming the source function known (cf. 
Mihalas 1978). In our work this is done by a short characteris- 
tics method (Olson & Kunasz 1987). Irradiation of the disc by 
an external source with a given spectrum is accounted for by 
an appropriate boundary condition. 

2.2. Structure Equations 

In order to obtain the atomic population numbers in NLTE all 
processes populating or de-populating an atomic level are con- 
sidered. These are ionisation, recombination, excitation and de- 
excitation, caused either by radiation or collision. Each level / 
of each ion has one rate equation describing the modification 
of the population density «, with time f; 



(8) 
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Pij denotes the rate coefficients, consisting of radiative and col- 
lisional components. Since we assume a hydrostatic stratifica- 
tion, we have 



drii 
'dt 



= 0. 



(9) 



Assuming that the radial component of the gravitation of 
the central object equals the centrifugal force of the Keplerian 
rotation of the disc, the hydrostatic equilibrium of the disc only 
is determined by the vertical component of the gravitation 



(10) 



with P denoting the total (gas and radiation) pressure. 

Another fundamental equation of the vertical structure de- 
scribes the local energy balance: 



^mech 



(11) 



The viscously generated energy fimech is equal to the radiative 
energy loss the convective energy ficonv term is neglected 
in the following. For standard c-discs, the viscously generated 
energy reads 



and 



^rad = 47r J" (/7(v,z) - xiv,z)J(y,z))dv. 



Introducing now the mass column depth m as 

oo 

m{z) - J pdz 



(12) 



(13) 



(14) 



with the total column mass Mq at the midplane the energy dis- 
sipated at each depth is 



Emechim) = - —^v(m)p. 



(15) 



Here v(m) denotes the depth-dependent kinematic viscosity 

with f>0 (16) 



v(m) = av(f+l)(^j 

with the damping factor which has been introduced by Kriz 
& Hubeny (1986) to avoid numerical instabilities at the disc 
surface, v is the depth-averaged kinematic viscosity with 



Mo 



v(w) dm . 

V corresponds to v in Eq. (|3} and can be determined by 



V - 



Re 



Re 



(17) 



(18) 



with V0 denoting the Keplerian angular velocity and Re the ef- 
fective Reynolds number (Lynden-Bell & Pringle 1974). The 



prescription of the viscosity using the Reynolds number has 
the advantage that no further assumptions concerning first val- 
ues of disc height and speed of sound have to be made, as would 
be necessary using the c-prescription. Furthermore, it is possi- 
ble to describe the depth dependency of the viscosity using the 
column mass as depth variable. The solution of the energy bal- 
ance is obtained with a generalised Unsold-Lucy method (Lucy 
1964, Dreizler 2003). 

Finally, the total particle density consists of the sum of 
the population numbers n of the NLTE and LTE levels I of all 
ions / of all elements x plus the electron density 11^: 



element ion /NLTE LTE ^ 



+ rip 



x=\ (=1 V /=! 1=1 
The equation of charge conservation reads 



element ion 



/'NLTE LTE 
dl + 



V 1=1 



1=1 



(19) 



(20) 



with charge q(i) of the ion /. 

The solution of the system of equations consisting of the 
radiation transfer equation, the equations of energy balance 
and hydrostatic equilibrium, the rate equations and the equa- 
tions of charge and particle conservation is done in an itera- 
tive scheme, the so-called Accelerated Lambda Iteration (ALI, 
Werner & Husfeld 1985; Werner et al. 2003). The input param- 
eters of our models are mass and radius of the central object, 
radius of the disc ring, mass accretion rate and Reynolds num- 
ber. Furthermore, the atomic data and an appropriate frequency 
grid are specified (see e.g. Rauch & Deetjen 2003). 

2.3. LTE Start Models 

To avoid numerical instabilities at the beginning of the NLTE 
calculations it is necessary to establish suitable start models un- 
der the assumption of LTE. In the following, we summarise the 
main steps creating such models, according to Hubeny (1990). 

To determine the atomic LTE population numbers n the 
Boltzmann equation 



n: 



Li g-(E,-Ej)lkT 



V 8j 
and the Saha equation 



n 



up 



nio 



2 / ImnekT 
n 



■S"P „-(£u„-£iow)/<rr 



(21) 



(22) 



are used instead of the NLTE rate equations. Here, g denotes 
the statistical weights, E the excitation or ionisation energy, 
me the electron mass and h the Planck constant. Furthermore, 
in the case of LTE the source function Sy equals the Planck 
function By 



S y = By 



(23) 



To get an analytical expression for the vertical temperature 
structure one combines the first momentum of the specific in- 
tensity and the equation of the energy balance. Then the equa- 
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tion for the vertical temperature structure in the case of LTE 
finally reads 



i=0° 



r4 3. 



with £ - Kb Mq/t and 



Kb 



OO 

4> 



/p) By dv . 



(24) 



(25) 



The equation of the hydrostatic equilibrium has the same 
form as in the case of NLTE, but can alternatively be trans- 
formed into a differential equation of second order: 



d^P _ clGMi, 



(26) 



The solution of this equation is done numerically. The upper 
boundary condition reads, following Hubeny (1990), 



P(l) 



with 



miCj 1 



fix) = ^e^'kix) 



and 



k{x) 



CO 

- \ e '^dt. 
n J 



(27) 



(28) 



(29) 



Here, Hg and denote the gas and radiation pressure scale 
height, defined as 



Icl 



0- 4 GU^ 



(30) 



(31) 



with Cg denoting the sound speed associated with gas pressure 

(32) 



cl = ^ 



2.4. The Synthetic Spectrum 

Having now calculated the vertical structures and spectra of the 
individual disc rings, the ring spectra are integrated to get the 
spectrum of the whole accretion disc: 



So 2k 



I{v, i) - cos(/) 



r) rd(f>dr . 



(33) 



Ri 



Here, Ri and Ro denote the inner and outer radius of the disc, 
; is the inclination angle (cf. Fig. O. The integration over the 
azimuthal angle <p is done in intervals of 1°. In this last step, 
spectral lines become broadened due to the Keplerian rotation 
of the disc. 



inclination angle 



Observer 




i=90 



central object 



accretion disk 
Fig. 3. Accretion disc geometry. 

3. Synthetic Spectrum for AMCVn 

In order to test and gain experience with AcDc our first applica- 
tion was the calculation of the synthetic spectrum of AM CVn 
for a comparison with results of a recent analysis performed by 
Nasser et al. (2001). 

3.1. AMCVn 

AMCVn is the prototype of the so-called AMCVn stars 
or helium cataclysmics, a subgroup of the cataclysmic vari- 
ables. They are thought to be the end product of the evolu- 
tion of close binary systems (El-Khoury & Wickramasinghe, 
2000). The first observations of AMCVn have been made by 
Malmquist (1936) and Humason & Zwicky (1947). Greenstein 
& Matthews (1957) classified AMCVn as a helium-rich white 
dwarf, Burbidge et al. (1967) as a quasi-stellar object, and 
Wambler (1967) as a hot star After the discovery of periodic 
variability in the light curve (Smak 1967), Warner & Robinson 
(1972) proposed AM CVn to be a close binary system with on- 
going mass transfer. Right now, these systems are believed to 
be interacting white dwarf binary systems, consisting of a de- 
generate C/O white dwarf as primary and a semi-degenerate 
low-mass secondary, composed of almost pure helium. The 
secondary fills its Roche volume and loses mass via Roche- 
Lobe overflow onto the primary, generating an accretion disc 
around it. 

AMCVn itself has been analysed by Nasser et al. (2001) 
using TLUSDISC (Hubeny 1990). They showed that the sys- 
tem consists of a primary of about 1 . 1 M© and a secondary of 
about 0.09 Mq. The radius of the primary is 4600 km, and the 
mass accretion rate is about 3 ■ lO^^Mo/yr. For the calcula- 
tions, we assumed an inner radius of 1 .4 R* and an outer ra- 
dius of 15 R*, both values were varied to explore the influence 
of the disc size onto the spectrum. The Reynolds number was 
set to 15 000, the damping factor ( was chosen to be 0.001. 
Convective energy transport is not yet included in AcDc, so we 
had to neglect convection, fortunately without getting numer- 
ical problems in the outer part of the disc. The number ratio 
H/He was set to 10"^, and abundances of carbon, nitrogen, oxy- 
gen, and silicon were assumed to be solar. The self-consistent 
inclusion of metals is an improvement over the Nasser et al. 
(2001) models. Some details concerning the ions, levels and 
lines we used in our calculations are shown in Table 1. The 
atomic data are taken from the opacity project (Seaton et al. 
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Table 1. Some details concerning the ions and the number of 
NLTE levels and lines used in our calculations. 



Ion 


NLTE levels 


lines 


Ion 


NLTE levels 


lines 


Hi 


16 


29 


Nil 


2 





Hn 


1 


- 


Niii 


34 


67 


He I 


44 


31 


Niv 


34 


53 


Hen 


32 


59 


Nv 


36 


56 


He III 


1 




Nvi 


1 





Cii 


16 


25 


On 


26 


36 


Cm 


58 


124 


Oiii 


28 


37 


Civ 


36 


86 


Oiv 


11 


5 


Cv 


1 





Ov 


6 


1 








Ovi 


36 


48 








Ovii 


1 
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Fig. 4. Optical spectra of three accretion disc models with outer 
radii of 11 (soHd Hne), 13 R* (dotted) and 15 R* (dashed). 
The mass accretion rate is 3 ■ 10"^ Mo/yr and the incUnation is 
10°. 



1994) and the Kurucz line lists (1991). The Lyman, B aimer and 
Paschen series of H i, the Lyman series of He: and the Lyman, 
Balmer, Paschen and Bracket series of Hen, some Lyman and 
Balmer lines of metals as well as resonance lines are Stark 
broadened, all other line profiles are Doppler broadened. For 
the line broadening of H we use VCS tables (Lemke 1997), 
for the line broadening of He: we use BCS tables (Barnard, 
Cooper & Shamey 1969) and Griem tables (Griem 1974), for 
He II we use VCS tables (Schoning & Butler 1989). In total, we 
calculated 38 individual disc rings. 

3.2. Influence of Disc Parameters on the Spectrum 

First, we examined the influence of different inner and outer 
radii onto the spectrum of the disc. We varied the inner ra- 
dius from 1.4 to 2R* and the outer radius from 11 to 15R*. 
Figurel^shows the optical spectrum of three disc models with 
outer radii of 1 1 R*, 13 R* and 15 R*. The mass accretion rate 
is 3 ■ 10"^ Mo/yr, the inclination is 10°. One can clearly see the 
increasing total flux with increasing outer radius because of the 
increasing radiating surface. Figure [S] shows details of the nor- 
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Fig. 5. Detail of Fig. 0] with three accretion disc models with 
outer radii of IIR* (solid Une), 13 R* (dotted) and 15 R* 
(dashed). 
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Fig. 6. Optical spectra of an accretion disc model seen under 
three different inclination angles, 10° (solid line), 36° (dotted) 
and 60° (dashed). The inner radius of the disc is 1.4R*, the 
outer radius 15 R* and the mass accretion rate is 3 ■ lO"** Mo/yr. 



malised spectrum. The line cores of He i become deeper with 
increasing outer radius because the large outer regions of the 
disc are cool enough to show strong He i lines. Variation of the 
inner radius shows similar, but smaller effects on the lines. 

Second, we analysed the influence of different inclination 
angles on the spectrum. Figure|6lshows a disc model seen un- 
der three different inclination angles. The total flux decreases 
with increasing inclination angle because of the decreasing pro- 
jected surface. As shown in Fig. narrow lines at small in- 
clination angles become broad lines at high inclination angles 
because of the increasing rotational broadening. 

3.3. Vertical Structure 

In Fig. Othe vertical distribution of temperature, density and 
Rosseland optical depth of three representative disc rings at 
1.4 R*, 7R* and 15 R* is shown. According to Eq. (0) disc 
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4700 



Fig. 7. Detail of Fig.|6lwith an accretion disc model seen under 
three different inclination angles, 10° (solid line), 36° (dotted) 
and 60° (dashed). 




-2 -1 
log m [g/cm^] 

Fig. 8. Vertical distribution of temperature, density and 
Rosseland optical depth of three representative disc rings with 
radii of 1 .4 R* (soHd line), 7 R* (dotted hne) and 15 R* (dashed 
line). 



rings from the outer part of the disc are cooler than inner disc 
rings. The temperature decreases monotonously from the mid- 
plane to the surface of the rings. Disc rings composed only of 
H and He, but without C, N, and O show an increase of the 
temperature in the outer layers near the surface. 



Fig. 9. Comparison of the observed spectrum (thin line) with 
model spectra of the AM CVn disc for 3 different outer radii. 




inclination =70°, 48°, 23° 



4000 4200 4400 4600 4800 5000 

Fig. 10. Comparison of the observed spectrum (thin line) with 
model spectra of the AM CVn disc for 3 different inclination 
angles. 

3.4. Model vs. Observation 

Finally we compared our accretion disc model spectra with an 
observed spectrum of AM CVn, obtained at the 6 m BAT of the 
SAO in May 2002. 

Figure |9] shows the observed spectrum and three model 
spectra with outer radii of 1 1 R* , 1 3 R* and 1 5 R* , respectively. 
The inner radius is 1.4R*, the inclination is 36° and the mass 
accretion rate is 3 ■ 10"'' Mo/yr. Owing to the larger and at the 
same time cooler radiating surface, a larger outer radius leads 
to an increase of the spectral line strengths of neutral helium 
compared to those of ionized helium. 

Figure [TUl shows the influence of the inclination angle on 
the spectrum. The disc models with inclination angles of 23°, 
48° and 70° extend from 1.4R* to 13 R* and the mass accre- 
tion rate is 3 ■ 10"^ Mo/yr. The larger the inclination angle, the 
stronger the spectral lines are broadened due to the increasing 
radial component of the Kepler velocity. 

In Fig. [^our best fit is shown. The disc model extends 
from 1.4R* to 13 R*, the inclination is 36° and the mass accre- 
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4000 4200 4400 4600 4800 500 

Fig. 11. Comparison of the observed spectrum (thin Hne) with 
a model spectrum of the AM CVn disc. 



tion rate is 3 ■ 10 'Mo/yr. The shapes of the Hei lines are in 
good agreement with the observation. 

4. Conclusions 

1 . We developed the new code AcDc for detailed NLTE cal- 
culations of accretion disc spectra of cataclysmic variables 
and compact X-ray binaries, solving the radiation transfer 
equation self-consistently together with the structure equa- 
tions under consideration of full metal Une blanketing as 
well as irradiation of the accretion disc by the central ob- 
ject. 

2. Variation of the extension of the disc has a significant in- 
fluence on the spectrum of the AM CVn disc. Lines of He i 
become deeper with increasing outer radius. In addition, 
the inclination angle of the disc has a clear influence on the 
spectral lines by Kepler rotation. 

3. The comparison of model spectra with an observed spec- 
trum of AM CVn yields 1.4R* for the inner radius and 
13 R* for the outer radius, the inclination is found to be 
36° and the mass accretion rate 3 ■ lO^^Mo/yr. This is in 
agreement with results of Nasser et al. (2001). 

For future work, we plan to include further physical processes 
in the program, e.g. Comptonisation as well as convective en- 
ergy transport, into the program. Furthermore, we have started 
to study winds from the accretion disc. A depth dependent 
atomic model will be employed to overcome numerical insta- 
bilities in irradiated discs due to extremely weak populated 
atomic levels at some depths. 
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